/** \file modelos.cpp
 *  \brief Implementação das funçoes dos modelos de diversos sistemas dinâmicos estudados durante a elaboração deste trabalho
 *  \author Geraldo Francisco de S. Rebouças | gfs.reboucas@gmail.com
 *  \since 28/11/2012
 *  \version 1.0
 *  \date 02/12/2012
 *  \see Veja a bibliografia deste trabalho */

#include "modelos.hpp"

vetor Duffing(double T, vetor X, double PV)
{/** \param T Recebe a variável independente, normalmente o tempo
     \param X Recebe o vetor de estados
     \param PV Recebe o parâmetro que desejamos variar durante as simulações, muito útil quando traçamos o diagrama de bifurcação do sistema    */

    if(X.retDim() != 2)
    {
        cerr << "Duffing \nDimensoes incompativeis \n";
        abort();
    }

    vetor xp(2, false);

    double inc = 0.1; // Incerteza dos parâmetros, 0.4 é legal
    double w = Rand(-inc, inc);

    double a1 = +0.2*(1.0 + w); // constante de amortecimento
    double a2 = -1.0*(1.0 + w); // constantes da mola não-linear
    double a3 = +2.0*(1.0 + w); // idem
    double a4 = +0.2*(1.0 + w); // Amplitude do forçamento
    double a5 = +1.0; // frequencia de excitação senoidal
    double u = a4*cos(a5*T); // Entrada do sistema

    xp(0) = X(1); // Sistema escrito com variáveis de estado
    xp(1) = -a1*X(1) - a2*X(0) - a3*pow(X(0), 3.0) + u;

    return xp;
}

vetor RCSJJ(double T, vetor X, double PV)
{/** \param T Recebe a variável independente, normalmente o tempo
     \param X Recebe o vetor de estados
     \param PV Recebe o parâmetro que desejamos variar durante as simulações, muito útil quando traçamos o diagrama de bifurcação do sistema    */

    if(X.retDim() != 2)
    {
        cerr << "RCSJJ \nDimensoes incompativeis \n";
        abort();
    }

    vetor xp(2, false);

    double inc = 0.1; // Incerteza dos parâmetros
    double w = Rand(-inc, inc);
    double b, a1, a2, a3; // Constantes do sistema
    b  = 1.0*(1.0 + w); // Valor nominal mais uma incerteza entre (-var,var)
    a1 = 1.0*(1.0 + w);
    a2 = 1.0*(1.0 + w);
    a3 = 0.5; // frequencia de excitação, \f$ \omega \f$

    double u = PV + 1.6*sin(T*a3); // Entrada do sistema

    xp(0) = X(1); // Sistema escrito com variáveis de estado
    xp(1) = b*u - a1*sin(X(0)) - a2*X(1);

    return xp;
}

vetor ROV(double T, vetor X, double PV)
{
    if(X.retDim() != 2)
    {
        cerr << "ROV \nDimensoes incompativeis \n";
        abort();
    }

    double inc = 0.25; // Incerteza dos parâmetros
    double w = Rand(-inc, inc);
    double m, C, A, rho; // Constantes do sistema
    m = 50.0*(1.0 + w); // Valor nominal mais uma incerteza entre (-var,var)
    C = 2.0*(1.0 + w);
    A = 0.25;
    rho = 1000.0;

    vetor xp(2, false);

    xp(0) = X(1); // Sistema escrito com variáveis de estado
    xp(1) = (PV - 0.5*rho*C*A*X(1)*fabs(X(1)))/m;

    return xp;
}

vetor RLD(double T, vetor X, double PV)
{
    if(X.retDim() != 2)
    {
        cerr << "RLD \nDimensoes incompativeis \n";
        abort();
    }

    double inc = 0.5; // Incerteza dos parâmetros
    double g = Rand(-inc, inc);
    double L, E, w, R, E0, Cj, Cd, a, b;
    L = 100.0e-6*(1.0 + g);
    E = PV; //80.0e-3; //2.0;
    w = 1400.0*M_PI*1.0e3;
    R = 60.0*(1.0 + g);
    E0 = 0.1*(1.0 + g);
    Cj = 400.0e-12*(1.0 + g);
    Cd = 0.1e-6*(1.0 + g);
    a = (Cj - Cd)/(2*Cd*Cj);
    b = (Cj + Cd)/(2*Cd*Cj);

    vetor xp(2, false);

    xp(0) = X(1);
    xp(1) = (E*sin(w*T) - R*X(1) - a*fabs(X(0)) - b*X(0) - E0)/L;

    return xp;
}

vetor VdP(double T, vetor X, double PV)
{
    if(X.retDim() != 2)
    {
        cerr << "VdP \nDimensoes incompativeis \n";
        abort();
    }

    double inc = 0.5; // Incerteza dos parâmetros
    double w = Rand(-inc, inc);
    double mu, A, B, C, f;
    mu = 1.0*(1.0 + w);
    A = 1.0*(1.0 + w);
    B = 1.0*(1.0 + w);
    C = 0.0*(1.0 + w);
    f = PV;

    vetor xp(2, false);

    xp(0) = X(1);
    xp(1) = mu*(1.0 - A*X(0)*X(0))*X(1) - B*X(0) - C + f;

    return xp;
}

vetor DuffingEst(double T, vetor X, double PV)
{/** \param T Recebe a variável independente, normalmente o tempo
     \param X Recebe o vetor de estados
     \param PV Recebe o parâmetro que desejamos variar durante as simulações, muito útil quando traçamos o diagrama de bifurcação do sistema    */

    if(X.retDim() != 2)
    {
        cerr << "DuffingEst \nDimensoes incompativeis \n";
        abort();
    }

    vetor xp(2, false);

    double a1 = +0.1; // constante de amortecimento
    double a2 = -1.0; // constantes da mola não-linear
    double a3 = +1.0; // idem
    double a4 = +2.0; // Amplitude do forçamento
    double a5 = +1.0; // frequencia de excitação senoidal
//    double u = a4*cos(a5*T); // Entrada do sistema
    double u = a4*sin(a5*T); // Entrada do sistema

// x'' + a1*x' + a2*x + a3*x^3 + a4*sin(a5*t)
    xp(0) = X(1); // Sistema escrito com variáveis de estado
    xp(1) = -a1*X(1) - a2*X(0) - a3*pow(X(0), 3.0) + u;

    return xp;
}

vetor RCSJJest(double T, vetor X, double PV)
{
    if(X.retDim() != 2)
    {
        cerr << "RCSJJEst \nDimensoes incompativeis \n";
        abort();
    }

    double b, a1, a2, a3;
    b  = 1.0; // Valor nominal das constantes
    a1 = 1.0;
    a2 = 1.0;
    a3 = 0.5; // Frequencia de excitação, \f$ \omega \f$

    double u = PV*sin(a3*T); // Entrada do sistema

    vetor xp(2, false);

    xp(0) = X(1);
    xp(1) = b*u - a1*sin(X(0)) - a2*X(1);

    return xp;
}

vetor RLDest(double T, vetor X, double PV)
{
    if(X.retDim() != 2)
    {
        cerr << "RLDEst \nDimensoes incompativeis \n";
        abort();
    }

    double L, E, w, R, E0, Cj, Cd, a, b;
    L = 100.0e-6;
    E = PV; //80.0e-3; //2.0;
    w = 1400.0*M_PI*1.0e3;
    R = 60.0;
    E0 = 0.1;
    Cj = 400.0e-12;
    Cd = 0.1e-6;
    a = (Cj - Cd)/(2*Cd*Cj);
    b = (Cj + Cd)/(2*Cd*Cj);

    vetor xp(2, false);

    xp(0) = X(1);
    xp(1) = (E*sin(w*T) - R*X(1) - a*fabs(X(0)) - b*X(0) - E0)/L;

    return xp;
}

vetor VdPest(double T, vetor X, double PV)
{
    if(X.retDim() != 2)
    {
        cerr << "VdPEst \nDimensoes incompativeis \n";
        abort();
    }

    double mu, A, B, C, f;
    mu = 1.0;
    A = 1.0;
    B = 1.0;
    C = 0.0;
    f = PV;

    vetor xp(2, false);

    xp(0) = X(1);
    xp(1) = mu*(1.0 - A*X(0)*X(0))*X(1) - B*X(0) - C + f;

    return xp;
}

vetor Lorenz(double T, vetor X, double PV)
{
    if(X.retDim() != 3)
    {
        cerr << "Lorenz \nDimensoes incompativeis \n";
        abort();
    }

    double a1, a2, a3;
    a1 = 16.0;
    a2 = 45.92;
    a3 = 4.0;

    vetor xp(3, false);
    xp(0) = a1*(X(1) - X(0));
    xp(1) =(a2 - X(2))*X(0) - X(1);
    xp(2) = X(0)*X(1) - a3*X(2);

    return xp;
}

matriz LorenzJac(double T, vetor X, double PV)
{
    if(X.retDim() != 3)
    {
        cerr << "LorenzJac \nDimensoes incompativeis \n";
        abort();
    }

    double a1, a2, a3;
    a1 = 16.0;
    a2 = 45.92;
    a3 = 4.0;

    matriz J(3, 3);

    J(0,0) = -a1;         J(0,1) = a1;      J(0,2) = 0.0;
    J(1,0) = a2 - X(2);   J(1,1) = -1.0;    J(1,2) = -X(0);
    J(2,0) = X(1);        J(2,1) = X(0);    J(2,2) = -a3;

    return J;
}

matriz LorenzVE(double T, vetor X, matriz Fi, double PV)
{
    if((X.retDim() != 3) or (Fi.retc() != 3) or (Fi.retl() != 3))
    {
        cerr << "LorenzVE \nDimensoes incompativeis \n";
        abort();
    }

    matriz Fip(3,3);
    Fip = LorenzJac(T, X, PV)*Fi;

    return Fip;
}
